1. Basic deterministic under-sampling

In this first notebook, we start with simple deterministic under-sampling using either a low-frequency or a high-frequency mask. Hence we want to generate a sampling mask $\Omega\in\{0,1\}^N$ such that $y = F_{\Omega} x$ where $x\in R^{N}$ is the original image. The data $y$ is then corrupted by additive noise $b$:

$$ y = F_{\Omega} x + b $$

where $b$ is a complex-valued circular Gaussian noise of variance $\sigma^2$ and $F_{\Omega} = \Omega F$ in this Cartesian setting (i.e. sampling on the grid). We consider:

1. a low-frequency Cartesian mask $\Omega$ is defined by the central lines of k-space.
2. a low-frequency Cartesian mask $\Omega$ is defined by a square box in k-space centered at $(k_x,k_y)=(0,0)$.
3. a high-frequency Cartesian mask $\Omega$ is defined as the complementary set of the lines defined in to 1)
4. a high-frequency Cartesian mask $\Omega$ is defined as the complementary set of the box defined in 2)

Based on these masks we generate the corresponding measured data from a reference MR image and finally perform Cartesian image reconstruction $\widehat{x}$ from the under-sampled data $y$ as follows:

$$ \widehat{x} = \tilde{F}^* y $$

where $F^*$ is the zero-filled inverse Fourier transform (ifft).

• Author: Philippe Ciuciu (philippe.ciuciu@cea.fr)
• Date: 02/03/2025
• Target: IEEE EMBS-SPS Summer School on Novel acquisition and image reconstruction strategies in accelerated Magnetic Resonance Imaging

#DISPLAY BRAIN PHANTOM
'''
在医学成像（如 MRI、CT）和信号处理中，幻影（Phantom）是指用于测试、校准或研究的人工数据或物理模型。
它通常是一个合成（模拟）图像或实验装置，用来模拟真实世界的结构，比如大脑、心脏或其他人体组织。
'''
# 1. 导入必要的库

%matplotlib inline
import sys
import os.path as op
import os
import math
import cmath

import numpy as np

import matplotlib.pyplot as plt

from skimage import data, io     # 用于图像处理。
import brainweb_dl as bwdl       # 用于下载和处理脑部 MRI 幻影数据。

# 2. 设置 Matplotlib 显示参数
plt.rcParams["image.origin"]="lower"    # 设置图像的原点在左下角（默认是左上角）。
plt.rcParams["image.cmap"]='Greys_r'    # 设置图像颜色映射为灰度反转（黑白）。

# 3. 下载并处理 MRI 图像
mri_img = bwdl.get_mri(4, "T1")[70, ...].astype(np.float32)    
'''
bwdl.get_mri(4, "T1") : 从 BrainWeb 数据库下载一个 T1 加权 MRI 扫描
70 : 获取 MRI 体数据（3D）中的第 70 层切片，即 2D 切片图像。
.astype(np.float32) : 转换数据类型为 float32，以便进行计算。
''' 
#mri_img = bwdl.get_mri(4, "T2")[120, ...].astype(np.float32)

# 4. 显示 MRI 图像
print(mri_img.shape)            # (256, 256)
img_size = mri_img.shape[0]     # 获取图像的宽度（假设是正方形图像）。

# 5. 绘制 MRI 图像
plt.figure()
plt.imshow(abs(mri_img))        # 显示 MRI 图像，abs() 主要用于防止负像素值影响显示。
#plt.imshow(abs(mri_2D),origin="lower")     # lower : 图像的原点（0,0）位于左下角（默认）。
#plt.imshow(abs(mri_2D),origin="upper")     # upper : 图像的原点（0,0）位于左上角（常见于 Matplotlib 默认设置）
plt.title("Original brain image")
plt.show()

(256, 256)

• To get started, the objective is to build a low-frequency sampling mask that consists of the central k-space lines
• Fourier transform the reference MR image to compute k-space data, add zero-mean Gaussian complex-valued noise
• Mask the data using the above defined mask
• Next, perform zero-filled MR image reconstruction from masked k-space data
• Study the impact of the number of lines in the mask on the final image quality

# Generate trivial Cartesian sampling masks
#img_size = 512

# 1️. 生成低频 k-space 采样掩码

mask="low_res"
factor_ctr = 4

#Objective: construct a k-space mask that consists of the central lines
# 0 entries: points not sampled, 1 entries: points to be sampled
# Initialize k-space with 0 everywhere and then place the "1" appropriately
kspace_masklines = np.zeros((img_size,img_size), dtype="float64")   # 创建一个 k-space 采样掩码（全 0）

low_res_size = img_size // factor_ctr + 1   # 计算低频中心区域的大小
idx_vec = np.linspace(img_size // 2 - low_res_size // 2, img_size // 2 + low_res_size // 2, low_res_size)  # 选取中心的低频 k-space 线
idx_vec_ =  idx_vec.astype("int")
print(idx_vec_)     # idx_vec_ 选取了第 96 到 160 行。

#Use fancy indexing to select low frequency lines only
kspace_masklines[idx_vec_, ] = 1  # 将中心区域的 k-space 线设置为 1（表示被采样）

'''
- 目标：在 k-space 里，仅采样低频区域的中心 k-space 线，模拟 Cartesian 采样策略（类似于 MRI 采样）。 
- 作用：只保留 k-space 中的中心部分信息，丢弃高频信息，以观察重建效果的影响。
'''

if 0:
    plt.figure()
    plt.imshow(kspace_masklines, cmap='Greys_r')
    plt.title("Sampling mask")
    plt.show()


# 2️. 定义 Fourier 变换和逆变换
norm = "ortho"   # 使用正交归一化的傅立叶变换
#norm = None

def fft(x):
    return np.fft.fft2(x, norm=norm)  # 执行 2D Fourier 变换，将图像转换为 k-space（MRI 采样域）。

def ifft(x):
    return np.fft.ifft2(x, norm=norm) # 执行 逆 Fourier 变换，从 k-space 还原图像。

# 3️. 计算 k-space 数据并加噪

# Generate the kspace data: first Fourier transform the image
kspace_data = np.fft.fftshift(fft(mri_img))     # 对 MRI 图像执行 Fourier 变换，获取 k-space 数据
#add Gaussian complex-valued random noise
signoise = 10       # 添加高斯复值噪声（独立添加到实部和虚部）
#kspace_data += np.random.randn(*mri_img.shape) * signoise * (1+1j)
# Simulate independent noise realization on the real & imag parts
kspace_data += (np.random.randn(*mri_img.shape) + 1j*np.random.randn(*mri_img.shape)) * signoise

# 4️. 使用低频采样掩码进行 k-space 采样
# Mask data to perform subsampling
kspace_data *= kspace_masklines     # 仅保留 k-space 低频部分（中心区域）

# 5️. 进行 MRI 重建
# Zero order solution  
# 零填充（Zero-filled）MRI 反变换重建
image_rec0 = ifft(np.fft.ifftshift(kspace_data))            # get back to the convention
'''
由于我们丢弃了高频信息，重建图像将缺失细节，但保留了大部分结构信息。
'''

# 6️. 可视化结果
fig, axs = plt.subplots(2, 2, figsize=(10, 10))
# 显示原始 MRI 图像
axs[0, 0].imshow(mri_img)
axs[0, 0].set_title("True image")
# 显示低频采样掩码
axs[0, 1].imshow(kspace_masklines)
axs[0, 1].set_title("Low frequency sampling mask")
# 显示被噪声影响的 k-space 数据
axs[1, 0].imshow(np.abs(kspace_data), vmax=0.01*np.abs(kspace_data).max())
axs[1, 0].set_title("k-space noisy data")
# 显示重建的图像
axs[1, 1].imshow(np.abs(image_rec0))
axs[1, 1].set_title("Zero-order recon")
plt.show()

[ 96  97  98  99 100 101 102 103 104 105 106 107 108 109 110 111 112 113
 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131
 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149
 150 151 152 153 154 155 156 157 158 159 160]

Q1  

• Question: based on the previous example, try to construct a low-freqency sampling mask defined by a square box centered around $(k_x, k_y) = (0,0)$.
• Then, replicate the same steps:
  1. Generate noisy masked data
  2. Perform zero-filled MR image reconstruction
  3. Visualize results and study the impact of both the noise level and the mask size

# Generate trivial Cartesian sampling masks
#img_size = 512

mask="low_res"
factor_ctr = 4

'''
img_size = 256：

factor_ctr = 2 → half_size = (256 // 2 + 1) // 2 = 64
factor_ctr = 4 → half_size = (256 // 4 + 1) // 2 = 32
factor_ctr = 8 → half_size = (256 // 8 + 1) // 2 = 16
'''

#Objective: construct a k-space mask that consists of the central lines
# 0 entries: points not sampled, 1 entries: points to be sampled
# Initialize k-space with 0 everywhere and then place the "1" appropriately
kspace_masklines = np.zeros((img_size,img_size), dtype="float64")

low_res_size = img_size // factor_ctr + 1
idx_vec = np.linspace(img_size // 2 - low_res_size // 2, img_size // 2 + low_res_size // 2, low_res_size)
idx_vec_ =  idx_vec.astype("int")
#print(idx_vec_)

#Use fancy indexing to select low frequency lines only
#kspace_masklines[idx_vec_, ] = 1
kspace_masklines[idx_vec_, idx_vec_[0]:idx_vec_[-1]] = 1    #  **square box**

if 0:
    plt.figure()
    plt.imshow(kspace_masklines, cmap='Greys_r')
    plt.title("Sampling mask")
    plt.show()

norm = "ortho"
#norm = None

def fft(x):
    return np.fft.fft2(x, norm=norm)

def ifft(x):
    return np.fft.ifft2(x, norm=norm)

# Generate the kspace data: first Fourier transform the image
kspace_data = np.fft.fftshift(fft(mri_img))
#add Gaussian complex-valued random noise
signoise = 10
#kspace_data += np.random.randn(*mri_img.shape) * signoise * (1+1j)
# Simulate independent noise realization on the real & imag parts
kspace_data += (np.random.randn(*mri_img.shape) + 1j*np.random.randn(*mri_img.shape)) * signoise
# Mask data to perform subsampling
kspace_data *= kspace_masklines

# Zero order solution
image_rec0 = ifft(np.fft.ifftshift(kspace_data))            # get back to the convention

fig, axs = plt.subplots(2, 2, figsize=(10, 10))
axs[0, 0].imshow(mri_img)
axs[0, 0].set_title("True image")
axs[0, 1].imshow(kspace_masklines)
axs[0, 1].set_title("Low frequency sampling mask")
axs[1, 0].imshow(np.abs(kspace_data), vmax=0.01*np.abs(kspace_data).max())
axs[1, 0].set_title("k-space noisy data")
axs[1, 1].imshow(np.abs(image_rec0))
axs[1, 1].set_title("Zero-order recon")
plt.show()

Remark 1

To construct a low-freqency sampling mask defined by a square box centered around $(k_x, k_y) = (0,0)$, simply replace

kspace_masklines[idx_vec_, ] = 1

with

kspace_masklines[idx_vec_, idx_vec_[0]:idx_vec_[-1]] = 1

# Generate trivial Cartesian sampling masks

'''
vector_factor_ctr = [2, 4, 8, 16]
vector_signoise = [5, 10, 15, 20]
'''

def Reconstruction_v1 (mask, factor_ctr, signoise):

    norm = "ortho"
    #norm = None
    def fft(x):
        return np.fft.fft2(x, norm=norm)
    def ifft(x):
        return np.fft.ifft2(x, norm=norm)

    kspace_masklines = np.zeros((img_size,img_size), dtype="float64")
    low_res_size = img_size // factor_ctr + 1
    idx_vec = np.linspace(img_size // 2 - low_res_size // 2, img_size // 2 + low_res_size // 2, low_res_size)
    idx_vec_ =  idx_vec.astype("int")

    if mask == "low_res" :
        
        #Use fancy indexing to select low frequency lines only
        kspace_masklines[idx_vec_, idx_vec_[0]:idx_vec_[-1]] = 1    #  **square box**

        # Generate the kspace data: first Fourier transform the image
        kspace_data = np.fft.fftshift(fft(mri_img))
        
        # Simulate independent noise realization on the real & imag parts
        kspace_data += (np.random.randn(*mri_img.shape) + 1j*np.random.randn(*mri_img.shape)) * signoise
        # Mask data to perform subsampling
        kspace_data *= kspace_masklines

        # Zero order solution
        image_rec0 = ifft(np.fft.ifftshift(kspace_data))            # get back to the convention

        return image_rec0

$\leadsto$   Here we write a function for the mask == "low_res" case (a low-freqency sampling mask defined by a square box centered around $(k_x, k_y) = (0,0)$), which generates the reconstructed image based on the given factor_ctr and signoise.

mask="low_res"
vector_factor_ctr = [2, 4, 8, 16]
vector_signoise = [0, 10, 20, 100]

fig, axs = plt.subplots(4, 4, figsize=(12, 12))
for i in range(len(vector_factor_ctr)) :
    factor_ctr = vector_factor_ctr[i]

    for j in range(len(vector_signoise)) :
        
        signoise = vector_signoise[j]

        image_rec0 = Reconstruction_v1 (mask, factor_ctr, signoise)
        axs[i, j].imshow(np.abs(image_rec0))
        # axs[i].set_title("test")
        # axs[i, j].set_title(f"factor_ctr = {factor_ctr}, signoise = {signoise}") 
        axs[i, j].set_xlabel(f"signoise = {signoise}")
        axs[i, j].xaxis.set_label_position('top')  # 移动 xlabel 到顶部
        axs[i, j].set_ylabel(f"factor_ctr = {factor_ctr}")


# Add a main title (suptitle) to the figure
fig.suptitle(f"Reconstruction Results for Mask: {mask}", fontsize=16)

# Display the figure with adjusted layout
plt.tight_layout(rect=[0, 0, 1, 0.96])  # Adjust layout to fit suptitle
plt.show()

$\leadsto$   In the case where mask="low_res", we plot 16 reconstructed images with factor_ctr varying in vector_factor_ctr = [2, 4, 8, 16] and signoise varying in vector_signoise = [0, 10, 20, 100] :

• Since only low-frequency components are preserved while high-frequency information is discarded, the reconstructed images retain the overall structure but lose fine details.
• Since the img_size = 256, as factor_ctr increases from $2 \to 4 \to 8 \to 16$, the length of square box (sampling region) decreases from from $128 \to 64 \to 32 \to 16$ respectively. This means that fewer low-frequency k-space samples are retained, leading to progressively lower image quality in the reconstruction.
• Changes in the noise level signoise have little impact on the reconstructed image quality, as the low-frequency components dominate the overall structure.

Solution given by Github

#Objective: construct a k-space mask consisting of a central box
# Initialize k-space with 0 everywhere and then place the "1" appropriately
kspace_maskbox = np.zeros((img_size,img_size), dtype="float64")
# use fancy indexing along rows
kspace_maskbox[idx_vec_, ] = 1

list_img_size = np.arange(0., img_size).tolist()
filtered_center = [x for x in list_img_size if x not in idx_vec_]
array_idx_center = np.array(filtered_center)
array_idx_center_ = array_idx_center.astype("int")
# use fancy indexing along cols
kspace_maskbox[:, array_idx_center_] = 0

# Generate the kspace data: first Fourier transform the image
kspace_data = np.fft.fftshift(fft(mri_img))
#add Gaussian complex-valued random noise
signoise = 10
# Simulate independent noise realization on the real & imag parts
kspace_data += (np.random.randn(*mri_img.shape) + 1j * np.random.randn(*mri_img.shape)) * signoise
# Mask data to perform subsampling
kspace_data *= kspace_maskbox

fig, axs = plt.subplots(2, 2, figsize=(10, 10))
axs[0, 0].imshow(mri_img)
axs[0, 0].set_title("True image")
axs[0, 1].imshow(kspace_maskbox)
axs[0, 1].set_title("Low frequency sampling mask")
axs[1, 0].imshow(np.abs(kspace_data),  vmax=0.01*np.abs(kspace_data).max())
axs[1, 0].set_title("k-space noisy data")
axs[1, 1].imshow(np.abs(image_rec0))
axs[1, 1].set_title("Zero-order recon")
plt.show()

Q2  

• Next question: Construct a k-space sampling mask that consists of the high-frequency lines, just discard the central lines
• Study the impact of the number of lines removed

# Generate trivial Cartesian sampling masks
#img_size = 512

mask="low_res"
factor_ctr = 4

'''
img_size = 256：

factor_ctr = 2 → half_size = (256 // 2 + 1) // 2 = 64
factor_ctr = 4 → half_size = (256 // 4 + 1) // 2 = 32
factor_ctr = 8 → half_size = (256 // 8 + 1) // 2 = 16
'''

#Objective: construct a k-space mask that consists of the central lines
# 0 entries: points not sampled, 1 entries: points to be sampled
# Initialize k-space with 0 everywhere and then place the "1" appropriately
kspace_masklines = np.zeros((img_size,img_size), dtype="float64")

low_res_size = img_size // factor_ctr + 1
idx_vec = np.linspace(img_size // 2 - low_res_size // 2, img_size // 2 + low_res_size // 2, low_res_size)
idx_vec_ =  idx_vec.astype("int")
#print(idx_vec_)

#Use fancy indexing to select low frequency lines only
kspace_masklines[idx_vec_, ] = 1
# kspace_masklines[idx_vec_, idx_vec_[0]:idx_vec_[-1]] = 1    #  **square box**
kspace_masklines = 1 - kspace_masklines

if 0:
    plt.figure()
    plt.imshow(kspace_masklines, cmap='Greys_r')
    plt.title("Sampling mask")
    plt.show()

norm = "ortho"
#norm = None

def fft(x):
    return np.fft.fft2(x, norm=norm)

def ifft(x):
    return np.fft.ifft2(x, norm=norm)

# Generate the kspace data: first Fourier transform the image
kspace_data = np.fft.fftshift(fft(mri_img))
#add Gaussian complex-valued random noise
signoise = 10
#kspace_data += np.random.randn(*mri_img.shape) * signoise * (1+1j)
# Simulate independent noise realization on the real & imag parts
kspace_data += (np.random.randn(*mri_img.shape) + 1j*np.random.randn(*mri_img.shape)) * signoise
# Mask data to perform subsampling
kspace_data *= kspace_masklines

# Zero order solution
image_rec0 = ifft(np.fft.ifftshift(kspace_data))            # get back to the convention

fig, axs = plt.subplots(2, 2, figsize=(10, 10))
axs[0, 0].imshow(mri_img)
axs[0, 0].set_title("True image")
axs[0, 1].imshow(kspace_masklines)
axs[0, 1].set_title("Low frequency sampling mask")
axs[1, 0].imshow(np.abs(kspace_data), vmax=0.01*np.abs(kspace_data).max())
axs[1, 0].set_title("k-space noisy data")
axs[1, 1].imshow(np.abs(image_rec0))
axs[1, 1].set_title("Zero-order recon")
plt.show()

Remark 2

To construct a k-space sampling mask that consists of the high-frequency lines, simply add

kspace_masklines = 1 - kspace_masklines

after

kspace_masklines[idx_vec_, ] = 1

# Generate trivial Cartesian sampling masks

'''
vector_factor_ctr = [2, 4, 8, 16]
vector_signoise = [5, 10, 15, 20]
'''

def Reconstruction_v2 (mask, factor_ctr, signoise):

    norm = "ortho"
    #norm = None
    def fft(x):
        return np.fft.fft2(x, norm=norm)
    def ifft(x):
        return np.fft.ifft2(x, norm=norm)

    kspace_masklines = np.zeros((img_size,img_size), dtype="float64")
    low_res_size = img_size // factor_ctr + 1
    idx_vec = np.linspace(img_size // 2 - low_res_size // 2, img_size // 2 + low_res_size // 2, low_res_size)
    idx_vec_ =  idx_vec.astype("int")

    if mask == "low_res" :
        
        #Use fancy indexing to select low frequency lines only
        kspace_masklines[idx_vec_, idx_vec_[0]:idx_vec_[-1]] = 1    #  **square box**


    elif mask == "high_res" :
        
        #Use fancy indexing to select low frequency lines only
        kspace_masklines[idx_vec_, ] = 1
        kspace_masklines = 1 - kspace_masklines     #  **high-frequency lines**

    # Generate the kspace data: first Fourier transform the image
    kspace_data = np.fft.fftshift(fft(mri_img))
        
    # Simulate independent noise realization on the real & imag parts
    kspace_data += (np.random.randn(*mri_img.shape) + 1j*np.random.randn(*mri_img.shape)) * signoise
    # Mask data to perform subsampling
    kspace_data *= kspace_masklines

    # Zero order solution
    image_rec0 = ifft(np.fft.ifftshift(kspace_data))            # get back to the convention

    return image_rec0
    

$\leadsto$   Here we write a function not only for the mask == "low_res" case, but also adding the mask == "high_res" case according to the above code block (a k-space sampling mask consists of the high-frequency lines), which generates the reconstructed image based on the given factor_ctr and signoise.

mask = "high_res"
vector_factor_ctr = [2, 4, 8, 16]
signoise = 10

fig, axs = plt.subplots(4, 1, figsize=(12, 12))
for i in range(len(vector_factor_ctr)) :
    factor_ctr = vector_factor_ctr[i]

    image_rec0 = Reconstruction_v2 (mask, factor_ctr, signoise)
    axs[i].imshow(np.abs(image_rec0))
    # axs[i].set_title("test")
    # axs[i, j].set_title(f"factor_ctr = {factor_ctr}, signoise = {signoise}") 
    axs[i].set_xlabel(f"signoise = {signoise}")
    axs[i].xaxis.set_label_position('top')  # 移动 xlabel 到顶部
    axs[i].set_ylabel(f"factor_ctr = {factor_ctr}")


# Add a main title (suptitle) to the figure
fig.suptitle(f"Reconstruction Results for Mask: {mask}", fontsize=16)

# Display the figure with adjusted layout
plt.tight_layout(rect=[0, 0, 1, 0.96])  # Adjust layout to fit suptitle
plt.show()

$\leadsto$   In the case where mask="high_res", we construct a k-space sampling mask that consists of the high-frequency lines, just discard the central lines. We plot 4 reconstructed images with factor_ctr varying in vector_factor_ctr = [2, 4, 8, 16] and signoise = 10 to study the impact of the number of lines removed :

• Since only high-frequency components are preserved while low-frequency information is discarded, the reconstructed images retain fine details but lose the overall structure.
• Since the img_size = 256, as factor_ctr increases from $2 \to 4 \to 8 \to 16$, the length (number) of the center mask line (excluded from sampling region) decreases from $128 \to 64 \to 32 \to 16$ respectively. This means that more high-frequency k-space samples are retained, leading to progressively sharper but noisier reconstructed images.
• The impact of changes in the noise level signoise on the reconstructed image quality will be discussed in the next question.

Solution given by Github

# 0 entries: points not sampled, 1 entries: points to be sampled
# Initialize k-space with "1" everywhere and then place the "0" appropriately

kspace_masklines = np.ones((img_size,img_size), dtype="float64")
low_res_size = img_size // factor_ctr + 1
idx_vec_upper = np.linspace(img_size // 2 - low_res_size // 2, img_size // 2 + low_res_size // 2, low_res_size)
idx_vec_upper_ =  idx_vec.astype("int")

# set central k-space lines to 0: use fancy indexing along rows
kspace_masklines[idx_vec_upper_, ] = 0

# Generate the kspace data: first Fourier transform the image
kspace_data = np.fft.fftshift(fft(mri_img))
#add Gaussian complex-valued random noise
signoise = 10
#kspace_data += np.random.randn(*mri_img.shape) * signoise * (1+1j)
# Simulate independent noise realization on the real & imag parts
kspace_data += (np.random.randn(*mri_img.shape) + 1j * np.random.randn(*mri_img.shape)) * signoise
# Mask data to perform subsampling
kspace_data *= kspace_masklines

# Zero order solution
image_rec0 = ifft(np.fft.ifftshift(kspace_data))            # get back to the convention

fig, axs = plt.subplots(2, 2, figsize=(10, 10))
axs[0, 0].imshow(mri_img)
axs[0, 0].set_title("True image")
axs[0, 1].imshow(kspace_masklines)
axs[0, 1].set_title("High frequency sampling mask")
axs[1, 0].imshow(np.abs(kspace_data),  vmax=0.01*np.abs(kspace_data).max())
#axs[1].imshow(np.abs(np.fft.ifftshift(kspace_data)), cmap='Greys_r')
axs[1, 0].set_title("k-space noisy data")
axs[1, 1].imshow(np.abs(image_rec0))
axs[1, 1].set_title("Zero-order recon")
plt.show()

Q3  

Question: Based on the previous example, try to construct a k-space mask that consists of the removing the central box centered in $(𝑘_𝑥,𝑘_𝑦)=(0,0)$.

• Then, replicate the same steps:
  1. Generate noisy masked data
  2. Perform zero-filled MR image recon
  3. Visualize results and study the impact of both the noise level and the mask size

# Generate trivial Cartesian sampling masks
#img_size = 512

mask="high_res"
factor_ctr = 4

'''
img_size = 256：

factor_ctr = 2 → half_size = (256 // 2 + 1) // 2 = 64
factor_ctr = 4 → half_size = (256 // 4 + 1) // 2 = 32
factor_ctr = 8 → half_size = (256 // 8 + 1) // 2 = 16
'''

#Objective: construct a k-space mask that consists of the central lines
# 0 entries: points not sampled, 1 entries: points to be sampled
# Initialize k-space with 0 everywhere and then place the "1" appropriately
kspace_masklines = np.zeros((img_size,img_size), dtype="float64")

low_res_size = img_size // factor_ctr + 1
idx_vec = np.linspace(img_size // 2 - low_res_size // 2, img_size // 2 + low_res_size // 2, low_res_size)
idx_vec_ =  idx_vec.astype("int")
#print(idx_vec_)

#Use fancy indexing to select low frequency lines only
#kspace_masklines[idx_vec_, ] = 1
kspace_masklines[idx_vec_, idx_vec_[0]:idx_vec_[-1]] = 1    #  **square box**
kspace_masklines = 1 - kspace_masklines

if 0:
    plt.figure()
    plt.imshow(kspace_masklines, cmap='Greys_r')
    plt.title("Sampling mask")
    plt.show()

norm = "ortho"
#norm = None

def fft(x):
    return np.fft.fft2(x, norm=norm)

def ifft(x):
    return np.fft.ifft2(x, norm=norm)

# Generate the kspace data: first Fourier transform the image
kspace_data = np.fft.fftshift(fft(mri_img))
#add Gaussian complex-valued random noise
signoise = 10
#kspace_data += np.random.randn(*mri_img.shape) * signoise * (1+1j)
# Simulate independent noise realization on the real & imag parts
kspace_data += (np.random.randn(*mri_img.shape) + 1j*np.random.randn(*mri_img.shape)) * signoise
# Mask data to perform subsampling
kspace_data *= kspace_masklines

# Zero order solution
image_rec0 = ifft(np.fft.ifftshift(kspace_data))            # get back to the convention

fig, axs = plt.subplots(2, 2, figsize=(10, 10))
axs[0, 0].imshow(mri_img)
axs[0, 0].set_title("True image")
axs[0, 1].imshow(kspace_masklines)
# axs[0, 1].set_title("Low frequency sampling mask")
axs[0, 1].set_title("High frequency sampling mask")
# axs[1, 0].imshow(np.abs(kspace_data), vmax=0.01*np.abs(kspace_data).max())
axs[1, 0].imshow(np.abs(kspace_data), cmap='gray', vmax=1*np.abs(kspace_data).max())
axs[1, 0].set_title("k-space noisy data")
axs[1, 1].imshow(np.abs(image_rec0))
axs[1, 1].set_title("Zero-order recon")
plt.show()

Remark 3

To construct a k-space mask that consists of the removing the central box centered in $(𝑘_𝑥,𝑘_𝑦)=(0,0)$, simply replace

kspace_masklines[idx_vec_, ] = 1

with

kspace_masklines[idx_vec_, idx_vec_[0]:idx_vec_[-1]] = 1

Then add

kspace_masklines = 1 - kspace_masklines

after

kspace_masklines[idx_vec_, idx_vec_[0]:idx_vec_[-1]] = 1

# Generate trivial Cartesian sampling masks

'''
vector_factor_ctr = [2, 4, 8, 16]
vector_signoise = [5, 10, 15, 20]
'''

def Reconstruction_v3 (mask, factor_ctr, signoise):

    norm = "ortho"
    #norm = None
    def fft(x):
        return np.fft.fft2(x, norm=norm)
    def ifft(x):
        return np.fft.ifft2(x, norm=norm)

    kspace_masklines = np.zeros((img_size,img_size), dtype="float64")
    low_res_size = img_size // factor_ctr + 1
    idx_vec = np.linspace(img_size // 2 - low_res_size // 2, img_size // 2 + low_res_size // 2, low_res_size)
    idx_vec_ =  idx_vec.astype("int")

    if mask == "low_res" :
        
        #Use fancy indexing to select low frequency lines only
        kspace_masklines[idx_vec_, idx_vec_[0]:idx_vec_[-1]] = 1    #  **square box**


    elif mask == "high_res" :
        
        ##Use fancy indexing to select low frequency lines only
        #kspace_masklines[idx_vec_, ] = 1
        #kspace_masklines = 1 - kspace_masklines     #  **high-frequency lines**

        #Use fancy indexing to select low frequency lines only
        #kspace_masklines[idx_vec_, ] = 1
        kspace_masklines[idx_vec_, idx_vec_[0]:idx_vec_[-1]] = 1    #  **square box**
        kspace_masklines = 1 - kspace_masklines

    # Generate the kspace data: first Fourier transform the image
    kspace_data = np.fft.fftshift(fft(mri_img))
        
    # Simulate independent noise realization on the real & imag parts
    kspace_data += (np.random.randn(*mri_img.shape) + 1j*np.random.randn(*mri_img.shape)) * signoise
    # Mask data to perform subsampling
    kspace_data *= kspace_masklines

    # Zero order solution
    image_rec0 = ifft(np.fft.ifftshift(kspace_data))            # get back to the convention

    return image_rec0

$\leadsto$   Here we write a function not only for the mask == "low_res" case, but also adding the mask == "high_res" case according to the above code block (a k-space mask consists of the removing the central box centered in $(𝑘_𝑥,𝑘_𝑦)=(0,0)$), which generates the reconstructed image based on the given factor_ctr and signoise.

mask="high_res"
vector_factor_ctr = [2, 4, 8, 16]
vector_signoise = [0, 10, 20, 100]

fig, axs = plt.subplots(4, 4, figsize=(12, 12))
for i in range(len(vector_factor_ctr)) :
    factor_ctr = vector_factor_ctr[i]

    for j in range(len(vector_signoise)) :
        
        signoise = vector_signoise[j]

        image_rec0 = Reconstruction_v3 (mask, factor_ctr, signoise)
        axs[i, j].imshow(np.abs(image_rec0))
        # axs[i].set_title("test")
        # axs[i, j].set_title(f"factor_ctr = {factor_ctr}, signoise = {signoise}") 
        axs[i, j].set_xlabel(f"signoise = {signoise}")
        axs[i, j].xaxis.set_label_position('top')  # 移动 xlabel 到顶部
        axs[i, j].set_ylabel(f"factor_ctr = {factor_ctr}")


# Add a main title (suptitle) to the figure
fig.suptitle(f"Reconstruction Results for Mask: {mask}", fontsize=16)

# Display the figure with adjusted layout
plt.tight_layout(rect=[0, 0, 1, 0.96])  # Adjust layout to fit suptitle
plt.show()

$\leadsto$   In the case where mask="high_res", we construct a k-space mask that consists of the removing the central box centered in $(𝑘_𝑥,𝑘_𝑦)=(0,0)$. We plot 16 reconstructed images with factor_ctr varying in vector_factor_ctr = [2, 4, 8, 16] and signoise varying in vector_signoise = [0, 10, 20, 100] to study the impact of both the noise level and the mask size :

• Since only high-frequency components are preserved while low-frequency information is discarded, the reconstructed images retain fine details but lose the overall structure.
• Since the img_size = 256, as factor_ctr increases from $2 \to 4 \to 8 \to 16$, the length (number) of the center mask line (excluded from sampling region) decreases from $128 \to 64 \to 32 \to 16$ respectively. This means that more high-frequency k-space samples are retained, leading to progressively sharper but noisier reconstructed images. In particular, for factor_ctr = 8, we observe artifacts around the outer contour, though they do not significantly affect diagnosis.
• Changes in the noise level signoise do have impact on the reconstructed image quality this time, as the high-frequency components dominate the fine details : As the noise level signoise increases, the quality of the reconstructed image progressively declines.

Solution given by Github

# 0 entries: points not sampled, 1 entries: points to be sampled
# Initialize k-space with "1" everywhere and then place the "0" appropriately
kspace_maskbox = np.ones((img_size,img_size), dtype="float64")
# use fancy indexing along rows
kspace_maskbox[idx_vec_, ] = 0

list_img_size = np.arange(0., img_size).tolist()
filtered_center = [x for x in list_img_size if x not in idx_vec_]
array_idx_center = np.array(filtered_center)
array_idx_center_ = array_idx_center.astype("int")
# use fancy indexing along cols
kspace_maskbox[:, array_idx_center_] = 1

# Generate the kspace data: first Fourier transform the image
kspace_data = np.fft.fftshift(fft(mri_img))
#add Gaussian complex-valued random noise
signoise = 10
#kspace_data += np.random.randn(*mri_img.shape) * signoise * (1+1j)
# Simulate independent noise realization on the real & imag parts
kspace_data += (np.random.randn(*mri_img.shape) + 1j*np.random.randn(*mri_img.shape)) * signoise
# Mask data to perform subsampling
kspace_data *= kspace_maskbox

# Zero order solution
image_rec0 = ifft(np.fft.ifftshift(kspace_data))

fig, axs = plt.subplots(2, 2, figsize=(10, 10) )
axs[0,0].imshow(mri_img, cmap='Greys_r')
axs[0,0].set_title("True image")
axs[0,1].imshow(kspace_maskbox, cmap='Greys_r')
axs[0,1].set_title("High frequency box sampling mask")
axs[1,0].imshow(np.abs(kspace_data),  cmap='gray', vmax=1*np.abs(kspace_data).max())
#axs[1].imshow(np.abs(np.fft.ifftshift(kspace_data)), cmap='Greys_r')
axs[1,0].set_title("k-space noisy data")
axs[1,1].imshow(np.abs(image_rec0), cmap='Greys_r')
axs[1,1].set_title("Zero-order recon")
plt.show()

Q4  

Comment on changing the under-sampling mask (band vs square) on the final image quality ?

# Generate trivial Cartesian sampling masks

'''
vector_factor_ctr = [2, 4, 8, 16]
vector_signoise = [5, 10, 15, 20]
'''

def Reconstruction (mask, mode, factor_ctr, signoise):
    '''
    mask : element in vector_mask = ["low_res", "high_res"]
    mode : element in vector_mode = ["band", "square"]
    factor_ctr : element in vector_factor_ctr = [2, 4, 8, 16]
    signoise : element in vector_signoise = [5, 10, 15, 20]
    '''

    norm = "ortho"
    #norm = None
    def fft(x):
        return np.fft.fft2(x, norm=norm)
    def ifft(x):
        return np.fft.ifft2(x, norm=norm)

    kspace_masklines = np.zeros((img_size,img_size), dtype="float64")
    low_res_size = img_size // factor_ctr + 1
    idx_vec = np.linspace(img_size // 2 - low_res_size // 2, img_size // 2 + low_res_size // 2, low_res_size)
    idx_vec_ =  idx_vec.astype("int")

    if mask == "low_res" and mode == "band":

        ##Use fancy indexing to select low frequency lines only
        kspace_masklines[idx_vec_, ] = 1

    elif mask == "low_res" and mode == "square":
        
        #Use fancy indexing to select low frequency lines only
        kspace_masklines[idx_vec_, idx_vec_[0]:idx_vec_[-1]] = 1    #  **square box**

    elif mask == "high_res" and mode == "band":

        ##Use fancy indexing to select low frequency lines only
        kspace_masklines[idx_vec_, ] = 1
        kspace_masklines = 1 - kspace_masklines     #  **high-frequency lines**

    elif mask == "high_res" and mode == "square":
        
        ##Use fancy indexing to select low frequency lines only
        #kspace_masklines[idx_vec_, ] = 1
        #kspace_masklines = 1 - kspace_masklines     #  **high-frequency lines**

        #Use fancy indexing to select low frequency lines only
        #kspace_masklines[idx_vec_, ] = 1
        kspace_masklines[idx_vec_, idx_vec_[0]:idx_vec_[-1]] = 1    #  **square box**
        kspace_masklines = 1 - kspace_masklines

    else:
        raise ValueError("Error: Check the input variables !")

    # Generate the kspace data: first Fourier transform the image
    kspace_data = np.fft.fftshift(fft(mri_img))
        
    # Simulate independent noise realization on the real & imag parts
    kspace_data += (np.random.randn(*mri_img.shape) + 1j*np.random.randn(*mri_img.shape)) * signoise
    # Mask data to perform subsampling
    kspace_data *= kspace_masklines

    # Zero order solution
    image_rec0 = ifft(np.fft.ifftshift(kspace_data))            # get back to the convention

    return image_rec0

$\leadsto$   Here is our final version of the Reconstruction(mask, mode, factor_ctr, signoise) function. This function can be used to simulate reconstructed images for either low-frequency (mask = "low_res") or high-frequency (mask = "high_res") sampling and for two different shapes : band (mode = band) or square (mode = square).

vector_mask = ["low_res", "high_res"]
vector_mode = ["band", "square"]
vector_factor_ctr = [2, 4, 8, 16]
vector_signoise = [0, 10, 20, 100]

for p in vector_mask :
    for q in vector_mode :

        mask = p
        mode = q

        fig, axs = plt.subplots(4, 4, figsize=(12, 12))
        for i in range(len(vector_factor_ctr)) :
            factor_ctr = vector_factor_ctr[i]

            for j in range(len(vector_signoise)) :
                
                signoise = vector_signoise[j]

                image_rec0 = Reconstruction (mask, mode, factor_ctr, signoise)
                axs[i, j].imshow(np.abs(image_rec0))
                # axs[i].set_title("test")
                # axs[i, j].set_title(f"factor_ctr = {factor_ctr}, signoise = {signoise}") 
                axs[i, j].set_xlabel(f"signoise = {signoise}")
                axs[i, j].xaxis.set_label_position('top')  # 移动 xlabel 到顶部
                axs[i, j].set_ylabel(f"factor_ctr = {factor_ctr}")


        # Add a main title (suptitle) to the figure
        fig.suptitle(f"Reconstruction Results for Mask: {mask} + Mode: {mode}", fontsize=16)

        # Display the figure with adjusted layout
        plt.tight_layout(rect=[0, 0, 1, 0.96])  # Adjust layout to fit suptitle
        plt.show()

Comment :

Personally,

• In the case where mask = "low_res", mode = "band" performs better than mode = "square" :

  For example, when factor_ctr = 8, which means the center sample band comes to be sufficient, we can still clearly see the structure using the band mode.

  In contrast, the reconstructed image becomes blurry when using the square mode under the same conditions.

• In the case where mask = "high_res", mode = "band" performs worse than mode = "square" :

  For example, when factor_ctr = 8, which means the center unsampled band is sufficient, the reconstructed image contains artifacts inside, causing the loss of fine details in the band mode.

  In contrast, if mode = "square", the artifacts appear around the outer contour, while the internal details are well preserved under the same conditions.